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Abstract.  The  construction  of  a  polynomial  interpolant  to  data  given  at 
finite  pointsets  in  ntJ  (or,  most  generally,  to  data  specified  by  finitely  many  linear 
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Abstract.  The  construction  of  a  polynomial  interpolant  to  data  given  at 
finite  pointsets  in  IR1*  (or,  most  generally,  to  data  specified  by  finitely  many  linear 
functionals)  is  considered,  with  special  emphasis  on  the  linear  system  to  be  solved. 
Gauss  elimination  by  segments  (i.e.,  by  groups  of  columns  rather  than  by  columns) 
is  proposed  as  a  reasonable  means  for  obtaining  a  description  of  all  solutions  and 
for  seeking  out  solutions  with  ‘good’  properties.  A  particular  scheme,  due  to  Amos 
Ron  and  the  author,  for  choosing  a  particular  polynomial  interpolating  space  in 
dependence  on  the  given  data  points,  is  seen  to  be  singled  out  by  requirements  of 
degree-reduction,  scale-invariance,  and  a  certain  orthogonality  requirement.  The 
close  connection,  between  this  particular  construction  of  a  polynomial  interpolant 
and  the  construction  of  an  H-basis  for  the  ideal  of  all  polynomials  which  vanish  at 
the  given  data  points,  is  also  discussed. 


1.  Introduction 

Polynomial  interpolation  in  d  variables  has  long  been  studied  from  the  fol¬ 
lowing  point  of  view:  one  is  given  a  polynomial  space,  F,  and  seeks  to  characterize 
the  point  sets  0  C  IRrf  for  which  the  pair  (0,  F)  is  correct  in  the  sense  that  F 
contains,  for  each  g  defined  at  least  on  0,  exactly  one  element  which  agrees  with 
g  on  0.  Further,  the  space  F  is  either  the  space  Ilfc  of  all  polynomials  on  IRrf  of 
degree  <  fc,  or  is,  more  generally,  a  D-invariant  space  spanned  by  monomials.  See 
R.  A.  Lorentz  [9]  for  an  up-to-date  accounting  of  the  many  interesting  efforts  in 
this  direction. 

In  contrast,  de  Boor  and  Ron  [5]  starts  with  an  arbitrary  finite  point  set 
0  C  IRd  and  proposes  a  particular  choice,  denoted  lie,  from  among  the  many 
polynomial  spaces  F  for  which  (0,  F)  is  correct.  To  be  explicit  (an  explanation 
can  be  found  at  the  end  of  Section  7),  the  definition  of  11©  makes  use  of  the 
least  term  of  a  function  /,  which,  by  definition,  is  the  first,  or  lowest- degree, 
nontrivial  term  in  the  expansion  of  /  into  a  sum  /  =  /o  +  /i  +  h  +  . . .  in  which 
fj  is  a  homogeneous  polynomial  of  degree  j,  all  j.  With  this,  II©  is  defined  as  the 
linear  span 

11©  :=  span{/i  :  /  €  expe} 

of  least  terms  of  elements  of  the  linear  span  expe  of  exponentials  e«  exp(0x), 
9  6  ©.  Because  of  this,  we  have  come  to  call  the  resulting  interpolant 

P&9 
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to  a  given  function  g,  defined  at  least  on  0,  the  least  interpolant  to  g. 

The  connection  between  the  least  interpolant  and  Gauss  elimination  (applied 
to  the  appropriate  Vandermonde  matrix)  was  already  explored  in  de  Boor  [2]  and 
de  Boor  and  Ron  [7].  However,  this  earlier  discussion  left  the  (wrong)  impression 
that  a  particular  variant  of  Gauss  elimination  was  needed  in  order  to  obtain  the 
least  interpolant.  The  present  paper  corrects  this  impression,  in  the  process  of 
examining  the  question  of  what  might  single  out  II©  from  among  all  polynomial 
spaces  F  for  which  (0,  F)  is  correct.  To  be  sure,  de  Boor  and  Ron  [8)  already 
contains  such  a  discussion,  but  not  at  all  from  the  point  of  view  of  elimination. 
As  in  [8],  the  discussion  here  is  actually  carried  out  in  the  most  general  situation 
possible,  namely  of  interpolation  to  data  of  the  form  A g,  with  A  an  arbitrary  linear 
map  to  IRn  from  the  space  n  of  all  polynomials  on  HT*.  Since  such  ‘data  maps’  are 
of  finite  rank,  it  would  be  trivial  to  extend  them  to  functions  on  JRrf  other  than 
polynomials. 

In  the  discussion,  two  properties  of  II©  are  seen  to  play  a  particularly  impor¬ 
tant  role:  (i)  P©  is  degree-reducing,  i.e., 

deg P©<7  <  deg g,  VgGlI; 

and  (ii)  II©  is  scale-invariant,  i.e.,  for  any  r  >  0  and  any  p  €  II©,  also  p(r)  €  II©. 
The  latter  property  is  equivalent  to  the  fact  that  II©  stratifies,  i.e.,  II©  is  the 
direct  sum  II©  =  IIq  k  of  its  homogeneous  subspaces 

n©,fc  :=  IT©  n  n£ 

(with  U°k  the  space  of  all  homogeneous  polynomials  of  exact  degree  k  (including 
the  zero  polynomial)).  Also,  the  most  intriguing  property  of  II©  (see  de  Boor  and 
Ron  [8]),  namely  that 

ne=  H  kerPt(£)> 

p|e=o 

is  looked  at  anew.  Here,  pt  is  the  leading  term  of  the  polynomial  p,  i.e.,  the  last, 
or  highest-degree,  nontrivial  term  in  the  expansion  p  =  po+Pi+P2+-ofp  with 
pj  homogeneous  of  degree  j,  all  j. 

In  addition,  it  is  shown  how  to  make  use  of  the  construction  of  n©  by  elim¬ 
ination  to  construct  an  H-basis  for  the  polynomial  ideal  of  all  polynomials  which 
vanish  on  0.  This  responds  to  questions  raised  by  algebraists  during  talks  on  the 
algorithm  given  in  [2]  and  [7]. 

The  paper  highlights  Gauss  elimination  by  segments.  In  this  generalization 
of  Gauss  elimination,  the  matrix  A  to  be  factored  is  somehow  segmented,  i.e.,  is 
given  in  the  form 

A  =  [Ao,Ai,...], 

with  each  segment,  Aj,  comprising  zero  or  more  consecutive  columns  of  A.  (Here 
and  throughout,  we  use  MATLAB  notation.)  Correspondingly,  elimination  is  to 
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proceed,  not  column  by  column,  but  segment  by  segment.  In  the  case  of  polynomial 
interpolation,  the  segmentation  naturally  arises  by  grouping  monomials  of  the 
same  degree  together  and  because  there  is  no  natural  ordering  for  the  monomials 
of  the  same  degree.  Because  of  this  particular  application,  we  start  the  indexing 
of  the  segments  of  A  with  0. 

By  definition,  Gauss  elimination  by  segments  applied  to  A  =  [Ao,  Atl.. .] 
produces  a  factorization 


■  •  • .]  =  A/[I?o,  •  •  •] 

with  M  invertible  and 


R  [Ho,  i?i, . . .] 

a  segmented  row-echelon  form,  the  segmentation  corresponding  to  that  of 
A.  This  means  that  R  is  block  upper  triangular,  R  —  (/?tJ)  say  with  R3  — 
[i?o j ;  /?  i  j ; . . .]  having  exactly  as  many  columns  as  Aj  does,  with  each  diagonal 
block  Rj j  onto.  (A  casual  search  of  the  literature  failed  to  turn  up  earlier  occur¬ 
rences  of  this  concept  of  a  segmented  row  echelon  form.  However,  in  hindsight,  the 
notion  is  so  natural  that  it  is  bound  to  have  been  used  before.) 

Except  for  trivial  cases,  a  segmented  matrix  has  infinitely  many  such  factor¬ 
izations.  Nevertheless,  all  such  factorizations  have  certain  properties  in  common. 
For  our  purposes,  the  most  important  property  is  the  fact  that,  for  each  j,  the 
row-space  of  Rjj  depends  only  on  A  (and  the  segmentation).  In  particular,  it  is 
independent  of  the  details  of  the  elimination  process  which  led  to  R.  We  take  the 
trouble  to  prove  this  since,  as  we  also  show,  the  row-space  of  R:j  provides  all  the 
information  needed  to  construct  the  space  11©  7  =  II©  n  n°  of  homogeneous  poly¬ 
nomials  of  degree  j  in  n©.  We  also  identify  the  numerical  procedure  proposed  in 
[2]  and  [7]  as  a  particularly  stable  way  for  obtaining  such  a  segmented  row-echelon 
form. 

For  simplicity,  the  paper  deals  with  real-valued  functions  only.  All  the  results 
are  true  if  matrix  transpose  is  replaced  by  conjugate  transpose  and,  correspond¬ 
ingly,  certain  quantities  by  their  complex  conjugate. 

Finally,  the  paper  is  also  meant  to  illustrate  the  perhaps  debatable  point  that 
it  is  possible  to  write  about  polynomial  interpolation  in  several  variables  without 
covering  entire  pages  with  formulae. 

The  paper  had  its  start  in  a  discussion  I  had  with  Nira  Dyn,  in  August  of  1992, 
concerning  possible  alternatives  to  the  algorithm  in  [7],  and  subsequent  discussions, 
in  February  of  1993,  concerning  possible  alternatives  to  the  least  interpolant.  In 
particular,  the  use  of  column  maps  below  is  my  response  to  Nira  Dyn’s  very  direct 
and  useful  way  of  describing  the  least  interpolant. 
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2.  Polynomials 


A  polynomial  in  d  variables  is  customarily  written  in  the  form 

P(x)  =  ^xac(a),  (2-1) 

a 


with 

•■•x(d)"(<<),  x  =  (®(1) . i(d))  6  lRd 

and  the  sum  taken  over  the  set 

TL%  =  {a  €  lRd  :  a(j)  £  2Z+,  all  j} 

of  multiindices.  Further,  the  coefficient  sequence  c  has  finite  support. 

Since  it  is  important  to  distinguish  between  the  polynomial  p  and  its  value 
p{x)  at  the  point  x,  a  notation  for  the  power  map  x  xa  is  needed.  There  being 
no  standard  notation  around,  we’ll  use  for  it  ()Q.  Here  is  the  formal  definition: 

()°  :  IRd  ->  1R  :  x  -4  xa. 

In  these  terms,  the  polynomial  p  of  (2.1)  can  be  written 

a 

As  already  indicated  by  the  unusual  positioning  of  the  coefficients,  this  formula 
for  p  can  be  viewed  as  the  result  of  applying  the  ‘matrix’ 

A':=[()“:a€ZZd] 

with  ‘columns’  ()“  to  the  ‘vector’  c  in  the  standard  way,  namely  by  multiplying 
the  ‘column’  with  index  a  with  the  entry  of  c  with  index  a  and  summing.  Thus, 

p  =  A'c. 

This  notation  (taken  from  de  Boor  [3])  turns  out  to  be  very  convenient  in  what  is 
to  follow. 

Formally,  we  think  of  X  as  a  linear  map,  defined  on  the  space 

dom  X  :=  (ZZd  — ►  IR)o  :=  {c  :  2Zd  ->  IR  :  #  supp  c  <  oo} 

of  finitely  supported  real- valued  sequences  indexed  by  Z2d  ,  and  mapping  into  the 
space  lRd  ->  IR  of  real-valued  functions  of  d  arguments. 
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Note  that  X  is  1-1  and  its  range  is  II,  the  space  of  polynomials  in  d  inde- 
terminates,  interpreted  here  as  functions  on  IR^.  Any  linear  subspace  F  of  IT  of 
dimension  n  is  of  the  form 

F  =  ran  AW, 

with  W  a  1-1  linear  map  from  ]Rn  into  domX,  hence  V  :=  XW  (or,  in  more 
traditional  language,  the  sequence  of  columns  of  V)  is  a  basis  for  F.  To  be  sure, 
both  W  and  V  are  linear  maps  on  lRn,  hence  are,  in  the  language  of  [3],  column 
maps.  This  language  is  meant  to  stress  the  fact  that  any  linear  map  23,  from  n- 
dimensional  coordinate  space  Fn  to  some  linear  space  S  over  the  scalar-field  F, 
is  necessarily  of  the  form 

B  =  (6l,...,6nJ:Fn^5:c^^6;C(», 

J 

with 

bj  :  =  Bij 

the  jth  column  of  B  and 

i;  :=  (0,.;.,q,  1,0,...). 

j- 1  zeros 

We  denote  by 

the  number  of  columns  in  the  column  map  B.  If  the  target,  5,  of  the  column  map 
B  €  L(F",S)  is  itself  a  coordinate  space,  S  =  Fm  say,  hence  B  is  (naturally 
identified  with)  an  m  x  n-matrix,  then  the  columns  of  B  as  a  column  map  are, 
indeed,  its  columns  as  a  matrix. 

For  any  linear  map  A,  whether  a  column  map  or  not,  we  denote  by 

ran  A,  ker  A 

its  range  (or  set  of  values)  respectively  its  kernel  or  nullspace. 

3.  Interpolation 

As  discussed  in  the  Introduction,  we  are  interested  in  interpolation  at  a  given 
finite  point  set  0  C  JR*1,  of  cardinality 

n  :=  #0, 

say.  Define 

A  :  II  — »  lRn  :  g  •-*  5|©  ~  (g(9)  :  0  €  Q). 
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Then,  /  €  IT  interpolates  to  g  at  0  iff  /  solves  the  equation 

A?  =  A  g.  (3.1) 

We  Y  we  introduced  the  map  A  here,  not  only  for  notational  convenience,  but 
also  in  order  to  stress  the  fact  that  most  of  the  discussion  to  follow  is  valid  for 
any  onto  linear  map  A  :  II  -+  IRn.  We  call  any  linear  map  A  :  II  -4  IRn  a  data 
map  and  note  that  it  is  necessarily  of  the  form  A  :<?►->•  (A tg  :  i  =  1,.. . ,n)  for 
some  linear  functionals  A,  on  n,  which  we  call  the  rows  of  A,  for  obvious  reasons. 
Given  an  onto  data  map  A,  we  say  that  (A ,F)  is  correct  (or,  F  is  correct  for 
A)  if  (3.1)  has  exactly  one  solution  in  F  for  every  g  g  II.  In  that  case,  we  denote 
the  unique  interpolant  in  F  to  g  g  II  by 

Pf9- 

Here,  for  completeness,  are  certain  known  facts  about  interpolation. 

Lemma  3.2.  Let  F  be  an  n-dimensional  polynomial  space,  let  V  g  L(fR.n ,  n)  be 
a  basis  for  F,  and  let  A  €  L(U,  Kln)  be  onto.  Then,  the  following  are  equivalent: 

(i)  (A,  F)  is  correct. 

(ii)  A  |F  is  invertible. 

(Hi)  A  is  1-1  on  F. 

(iv)  A  V  is  invertible. 

Further,  if  any  one  of  these  conditions  holds,  then 

PF  =  P(AP)-*A. 

Proof.  Since  A  is  onto,  i.e.,  its  range  is  all  of  IR",  correctness  of  (A ,F) 
is  equivalent  to  having  the  linear  map  F  -4  IR"  :  /  4  A/  be  invertible,  and, 
since  both  domain  and  target  of  this  map  have  the  same  dimension,  n,  this  is 
equivalent  to  having  this  map  be  1-1.  Finally,  since  V  is  a  basis  for  F,  the  linear 
map  P|F  :  IRn  -4  F  :  a  >-4  Va  is  invertible.  Since  AP  =  A|f  P|F,  it  follows  that 
A  |f  is  invertible  iff  AP  is  invertible.  □ 

In  order  to  investigate  what  spaces  F  are  correct  for  A,  we  look  at  all  possible 
polynomial  interpolants  to  a  given  g,  i.e.,  at  all  possible  solutions  of  (3-1).  Since 
X  is  1-1,  the  solutions  of  (3.1)  are  in  1-1  correspondence  with  solutions  of  the 
equation 

AX ?  =  Ag. 

This  equation  is  a  system  of  n  linear  algebraic  equations,  albeit  in  infinitely  many 
unknowns,  namely  the  power  coefficients  c(a),  a  g  7Z±,  of  the  interpolant.  Never¬ 
theless,  we  may  apply  Gauss  elimination  to  determine  all  solutions. 

4.  Gauss  elimination 

Gauss  elimination  (with  row  pivoting)  produces  a  factorization 

AX  =  MR, 
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with  M  invertible,  and  R  ‘in  row-echelon  form’,  i.e.,  ‘right-triangular’  in  the  follow¬ 
ing  sense.  Each  row  other  than  the  first  is  either  entirely  zero,  or  else  its  left-most 
nonzero  entry  is  strictly  to  the  right  of  the  left-most  nonzero  entry  in  the  preceding 
row.  The  left-most  nonzero  entry  (if  any)  in  a  row  is  called  the  pivot  element 
for  that  row,  and  the  row  is  called  the  pivot  row  for  the  unknown  associated 
with  the  column  in  which  the  pivot  element  occurs.  An  unknown  is  called  bound 
or  free  depending  on  whether  or  not  it  has  a  pivot  row.  It  is  a  standard  result 
that  an  unknown  is  bound  if  and  only  if  its  column  is  not  in  the  linear  span  of  the 
columns  to  the  left  of  it. 

Since  AX  is  onto  (X  being  onto  II  and  A  being  onto),  every  row  of  R  will  be 
a  pivot  row.  Let  (qi  , . . . ,  an)  be  the  corresponding  sequence  of  indices  of  bound 
unknowns,  and  set  F  =  ran  V  with 

V:=  ()"“]. 

Then  A l7  is  invertible,  hence  (A,  F )  is  correct. 

This  construction  is  quite  arbitrary.  As  should  have  been  clear  from  the 
language  used,  Gauss  elimination  depends  crucially  on  the  ordering  of  the  columns 
of  A.Y,  i.e.,  on  the  ordering  of  the  columns  of  X.  On  the  other  hand,  for  d  >  1, 
there  is  no  natural  ordering  of  the  columns  of  X,  i.e.,  of  the  monomials.  The  best 
we  can  do,  offhand,  is  to  order  the  monomials  by  degree.  Here,  to  be  sure, 

deg  Ac  =  max{|a|  :=  ||a||i  =  ]Ta(i)  :  c(a)  ^  0}. 


Equivalently,  degp  =  min{A: :  p  €  11*},  with 


Ilfc  :=  ran  A<)t, 


where 

and 


A<fc  :=  [Xq,  ■ .  • ,  X/t] 


Assume  that  X  is  graded,  i.e.,  segmented  by  degree,  i.e., 

X  =  [X0,Xl,...]. 


This  has  the  following  happy  consequence. 

Result  4.1  ([7;  (2.5)]).  If  X  is  graded,  then  polynomial  interpolation  from  the 
space  F  spanned  by  the  monomials  corresponding  to  the  bound  columns  of  AX  is 
degree-reducing,  i.e., 

deg Pf9  <  degp  Vp  €  II. 
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Proof.  Recall,  as  we  did  earlier,  that  an  unknown  is  bound  if  and  only 
if  its  column  is  not  a  linear  combination  of  preceding  columns.  In  particular,  any 
column  is  in  the  linear  span  of  the  bound  columns  not  to  the  right  of  it.  This 
implies  that,  if  deg  g  =  k  and  therefore  A g  €  ran  AA'c^,  then  A g  is  necessarily 
in  the  linear  span  of  the  bound  columns  in  \X<k-  Since  /  :=  Ppg  is  the  unique 
element  in  F  with  A /  =  A g,  it  follows  that  Ppg  is  already  in 

Ft  :=Fnnfc. 

□ 

In  particular,  for  d  =  I  and  A  =  je,  we  recover  in  this  way  the  standard 
choice  for  F,  namely  II<n.  However,  for  d  >  I,  F  still  depends  on  the  ordering  of 
the  columns  within  each  X *  and  there  is  no  natural  ordering  for  them.  It  is  for 
this  reason  that  we  now  consider  Gauss  elimination  by  segments. 

5.  Gauss  elimination  by  segments 

To  recall  from  the  Introduction,  Gauss  elimination  by  segments  applied  to 
the  segmented  matrix  A  =  [.4o,  Ai , . . .]  produces  a  factorization 

[Ao,  Ai, . . .]  =  M[Rq,  Ri,  ■ . .) 

with  M  invertible  and 

R  :=  [Ro.Fi,---] 

a  segmented  row-echelon  form,  the  segmentation  corresponding  to  that  of  A.  This 
means  that  R  is  block  upper  triangular,  R  =  {RtJ),  say,  with 

'Ujj' 

Rj  —  [i?oj !  R\j  »•••]  —  tyjj !  R-jj  >0]  Rjj  , 

0 

#Rj  =  #Aj,  and  with  each  diagonal  block  Rjj  onto.  (Use  of  semicolon  as  in 
MATLAB.) 

While,  except  for  trivial  cases,  a  segmented  matrix  has  infinitely  many  such 
factorizations,  some  properties  of  such  a  factorization  depend  only  on  A  and  the 
particular  segmentation  used.  Here  are  the  basic  facts. 

Lemma  5.1.  Let  A,  B,  C,  M,  S,  Q,  T  be  matrices  with 

A  =  [B,C]  =  Af  ^  q  £], 

M  invertible,  #B  =  #S,  and  S  onto.  Then 
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(i)  The  transpose  S'  of  S  is  (i.e.,  the  rows  of  S  form)  a  basis  for  the  row-space 
of  B,  i.e.,  for  ran  B1 . 

(ii)  The  row-space  ofT,  i.e.,  the  space  ran  T' ,  depends  only  on  A  and  #T  = 
if  A  -  #B.  Explicitly, 

ran T1  =  (y  €  Ht*T  :  (0 ,y)  €  ranA'}. 


(Hi)  For  all  c  1  ran  T1 ,  there  exists  b  so  that  ( b,c )  €  ker.4. 

Proof.  (i):  Since  #B  =  #S,  we  have  B  —  MS,  with  A/  invertible,  hence 
ran  B1  =  ran  S' M1  =  ran  S'. 

r  s'  o  i 

(ii):  Since  M  is  invertible,  ran  A'  =  ran  q,  j,,\.  In  particular,  for  any 
#T-sequence  y  and  any  #S-sequence  x,  (x,y)  6  ran  .4'  if  and  only  if 

(x.y)  =  (S'zl,Q'zl+Tlz2) 


for  some  sequence  z  =  (zi ,  ;■>).  In  particular,  since  S  is  onto  and  therefore  S'  is 
1-1,  such  x  is  zero  iff  —  0.  In  other  words,  (0,  y)  €  ran  .4'  if  and  onlv  if =  T'z> 
for  some  z2. 

(iii):  If  c  J.  ranT'  then  c  €  ker7\  therefore,  for  any  #S-vector  b,  .4(6,  c)  = 
A/(S6  -f  Qc,0).  It  follows  that  ( b,c )  €  ker.4  if  and  only  if  Sb  =  -Qc  and,  since  S 
is  onto  by  assumption,  such  a  choice  for  b  is  always  possible.  □ 

Corollary  5.2.  Let 

AX  -  MR 

with  M  invertible  and  R  =  (R, _,)  a  segmented  row-echelon  form,  segmented  cor¬ 
responding  to  the  segmentation  X  =  [A'o,  Xi,. . .)  of  X.  Then, 

(i)  for  any  k,  R'kk  is  a  basis  for 

{y  €  lRri  :  (0,y)  €  ran(A.Y<*)'}, 

with  rk  :=  #Ri<k-  In  particular,  ran  R'kk  depends  only  on  A. 

(ii)  for  any  fc,  =  rankflfcic  =  rankA.Y<i;  -  rankAX^fc. 

Proof.  (i):  This  follows  from  (ii)  of  Lemma  5.1,  with  the  choices  B  := 
AX<k,  C  :=  AXk,  S  (Rij  :  C  <  i,j  <  k),  Q  :=  [Rok; Rk~i,k],  and  T  Rkk. 
(ii):  Since  M  is  invertible, 


rankAXot  =  rank  R<k, 

and  the  latter  number  equals  J2j<k  #R'jj  s'nce  ^<fc  =  (^»>  :  i  =  0,1,... ;j  = 
0, . .  • ,  k)  is  block  upper  triangular  with  each  diagonal  block  onto.  In  particular, 
each  Rjj  has  full  row  rank,  hence  #R'jj  =  rank  Rjj,  all  j.  □ 

We  briefly  discuss  the  construction  of  a  segmented  row-echelon  form  for  a 
given  segmented  matrix  A  =  [Ao,  Ai, . . .]. 
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It  follows  from  Lemma  5.1  that,  in  step  0  of  Gauss  elimination  by  segments 
applied  to  the  segmented  matrix  A  =  [.do.  A> , . . .],  one  determines,  by  some  means, 
some  basis  R^  for  ran  Ag.  While  this  could,  of  course,  be  done  by  Gauss  elimina¬ 
tion,  with  the  nontrivial  rows  of  the  resulting  row-ecr  Aon-form  for  A0  providing 
the  desired  basis,  i.e.,  providing  the  rows  for  the  matrix  Roo,  more  stable  proce¬ 
dures  come  to  mind.  For  example,  one  could  construct  R^o  as  an  orthogonal  basis, 
or  even  an  orthonormal  basis,  for  the  row-space  of  Ao-  The  favored  numerical  pro¬ 
cedure  for  this  would  be  to  construct  the  QR  factorization  for  A'0,  preferably  using 
Householder  reflections  and  column  interchanges,  and  use  the  resulting  orthonor¬ 
mal  basis  for  ranAg. 

Since  R^0  is  a  basis  for  ran  Ag,  there  is  a  unique  No  so  that  Aq  =  NqRoo-  In 
fact,  N0,  or  some  left  inverse  for  N0,  is  usually  found  during  the  construction  of 
Roo  In  any  case,  since  the  rows  of  Roo  are  a  basis  for  the  row-space  of  A0,  No  is 
necessarily  1-1,  hence  can  be  extended  to  some  invertible  matrix  A/o  =  [Ag, . . .]. 
This  gives  the  factorization 


A  =  [Aq,  A\ ,  A2,  •  •  ■]  =  A/q 


Roo 

0 


Q 

T\  T2 


Subsequent  steps  concentrate  on  the  segmented  matrix  T  :=  [Ti.T-j,  •  •  •]. 

This  discussion  of  Gauss  elimination  by  segments  was  only  given  in  order  to 
indicate  the  many  different  ways  available  for  constructing  a  factorization  A  =  MR 
of  a  segmented  matrix  A  =  [Ao,  Ai , . . .]  into  an  invertible  M  and  a  correspondingly 
segmented  row-echelon  form.  The  existence  of  such  a  factorization  was  never  in 
doubt,  as  any  factorization  A  =  MR  obtained  by  Gauss  elimination  with  partial 
pivoting  is  seen  to  be  of  the  desired  form  after  R  is  appropriately  blocked. 


6.  A  “natural”  choice  for  F 

We  have  already  seen  that  a  correct  F  can  be  so  chosen  that  Pp  is  degree- 
reducing.  The  following  result  characterizes  all  such  choices  F. 

Result  6.1  ([8]).  If  (A,  F)  is  correct,  then,  Pf  is  degree-reducing  if  and  only  if 
each  F/t  =  F  fl  rifc  is  as  large  as  possible. 

We  give  a  simple  proof  of  this  in  a  moment,  for  completeness.  But  first  we 
explore  what  limits  if  any  might  actually  be  imposed  on  the  dimension  of  Ft.  Since 
correctness  of  (A,  F)  implies  that  A  is  1-1  on  F,  A  must  be  1-1  on  every  F*,  hence, 
for  every  k, 

dimFfc  =  dim  A(Ffc)  <  dim  AX<k  =  rankAAf<*. 

Moreover,  equality  is  possible  here  for  every  A;  it  is  achieved  by  the  subspace  F 
spanned  by  the  monomials  corresponding  to  bound  columns  in  A[Ao,  X\ , . . .].  This 
proves  the  following 
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Corollary  6.2.  Assume  that  (A,  F)  is  correct.  Then,  Pf  is  degree-reducing  if  and 
only  if 

dim  Fk  =  rank  AX<k,  k  =  0,1,2, —  (6.3) 

In  other  words,  if  V  =  [V0,  Vj, . . .]  is  a  graded  basis  for  F,  i.e.,  V<k  is  a  basis 
for  Fk  for  every  k,  then  Pf  cannot  be  degree-reducing  unless 

#F<fc  =  rankAX<tl  *  =  0,1,2,....  (6.4) 

In  terms  of  any  segmented  row-echelon  form  R  =  (R,j)  for  A[A'o,  X\, . . .],  (6.4)  is, 
by  (ii)  of  Corollary  5.2,  equivalent  to 

#14  =  #R'kk,  k  =  0,1,2,.... 

Note  that  consideration  of  a  graded  basis  for  F  imposes  no  restriction  since 
such  a  basis  can  always  be  constructed  inductively,  starting  with  *  =  —  1 :  With  a 
basis  V<fc  for  Fk  already  in  hand,  a  requisite  Vk+i  is  obtained  by  completing  V<k 
to  a  basis  for  Fk+i .  For  the  proof  of  Result  6.1,  we  require  the  evident  fact  that 
we  are  free  to  include  in  14+ i  any  one  element  from  Fk+l\Fk. 

Proof  of  Result  6.1  If  Pf  is  degree- reducing  and  G  is  any  polynomial 
space  correct  for  A,  then  Pf  is  1-1  on  G  (since,  by  Lemma  3.2,  Pp  =  V(AV)-lA 
with  both  V  and  (AV)-1  1-1,  and  also  A  is  1-1  on  G ),  therefore,  for  any  k, 
dim  Gk  =  dim  Pp(Gk)  <  dim  Ft,  since  Pp(Gk)  C  Pf(IUi)  C  II*.  D  F  =  Ft  (using 
the  fact  that  Ff(flt)  C  IJk  since  Pf  is  degree-reducing). 

Conversely,  if  Pf  fails  to  be  degree-reducing,  then  there  exists  g  with  k  := 
deg <7  <  deg  PfQ  =:  *'•  Since  PfQ  €  Ft>\Ft«_i,  we  may  (by  the  last  sentence  before 
this  proof)  include  it  in  a  graded  basis  V  for  F.  Let  U  be  the  column  map  obtained 
from  V  by  replacing  Ppg  by  g.  Then  AV  =  A U,  hence  also  AU  is  invertible  and 
therefore  also  G  :=  ran  U  is  correct  for  A.  However,  dirnG*.  >  dim  Ft.  □ 

We  are  now  ready  to  discuss  use  of  an  appropriately  segmented  row-echelon 

form 

R  =  («,,)  =  A/"1  A[*o,  *i,H 

for  AX  in  checking  whether  a  polynomial  space  F  satisfying  (6.3)  is  correct  for 
interpolation  at  A.  If  V  is  any  basis  for  F,  then  (A ,F)  is  correct  iff  the  matrix  AV 
is  invertible.  With  V  =:  XW ,  this  is  equivalent  to  having  RW  invertible  (since  M 
is  invertible  by  assumption).  Now  assume  that  V  =  [Vb,  V, . . .]  is  a  graded  basis 
for  F.  Let 

Vk=:XWk,  k  =  0,1 . 

Then  V,  =  £,<;  XtWtJ,  all  i,j,  i.e.,  the  column  map  W,  from  IRn  into  dom  X, 
given  by  V  =:  XW ,  is  1-1  and  block  upper  triangular,  with  Wij  having  #Ar,  rows 
and  #V)  columns.  Since  (6.3)  holds,  #V)  =  #F'  j,  all  j.  Hence,  the  matrix  RW  is 
block  upper  triangular,  with  square  diagonal  blocks.  Therefore,  RW  is  invertible 
if  and  only  if  RkkWkk  is  invertible  for  each  k.  We  have  proved  the  following. 
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Proposition  6.5.  Let  MR  be  any  factorization  for  A[Xo,  Xi, . . .]  with  M  invert¬ 
ible  and  R  =  (Rij)  a  correspondingly  segmented  row-echelon  form.  Assume  that 
the  polynomial  space  F  satisfies  (6.3),  and  let  V  =  XW  be  a  graded  basis  for  F, 
hence  W  =  (Wtl)  is  block  upper  triangular,  with  Wl}  having  #X,  rows  and  #R'jj 
columns.  Then  (A,  F)  is  correct  if  and  only  if  RkkWkk  is  invertible  for  every  k. 

Note  that  this  condition  only  involves  the  diagonal  blocks  of  W .  This  suggests 
that  we  choose  all  off-diagonal  blocks  of  W  to  be  zero,  thus  keeping  the  structure 
of  F  simple.  This  is  equivalent  to  the  requirement  that 

F  =  0  XkWkk, 

k 

which,  in  turn,  is  equivalent  to  having  F  be  scale-invariant. 

Note  further  that  the  condition  of  the  Proposition  can  always  be  met,  since, 
by  construction,  each  Rkk  is  onto,  hence  has  right  inverses. 

If  Rkk  is  square,  then  invertibility  of  RkkWkk  implies  that  Wkk  is  invertible, 
hence  then 

ran  A'fclVfcjfc  =  ran  Xk  =  n°. 

In  other  words,  there  is  then  no  real  choice;  regardless  of  how  Wkk  is  chosen,  F 
contains  all  homogeneous  polynomials  of  degree  k. 

If  Rkk  is  not  square,  then  the  invertibility  of  RkkWkk  does  not  completely 
determine  the  kth  homogeneous  component  of  F.  How  is  one  to  choose  from  among 
the  infinitely  many  possibilities? 

The  choice 

Wkk  =  R'kk,  k  =  0,1,2,...  (6.6) 

suggests  itself,  as  it  guarantees  that  Rkk  Wkk  is  invertible  (since  Rkk  is  onto).  More¬ 
over,  if  R'f.,.  was  chosen  as  an  orthonormal  basis,  then  the  diagonal  block  RkkWkk 
becomes  an  identity  matrix,  thus  facilitating  the  calculation  (by  backsubstitution) 
of  the  coefficient  vector  c  of  the  resulting  interpolant  Xc  from  F  to  the  given  g. 
Indeed,  the  interpolant  is  of  the  form  XWa,  hence  a  solves  the  linear  system 


A  XW?  =  A  g. 


therefore  solves  the  equivalent  linear  system 

RW1  =  M~lAg.  (6.7) 

From  a,  the  coefficient  vector  c  in  the  power  form  Xc  for  the  interpolant  is  obtained 
as 

c  =  Wa. 

The  resulting  interpolation  scheme  Pp  has  additional  nice  properties.  How¬ 
ever,  we  don’t  bother  to  derive  them  here  since  Pf  has  a  major  flaw:  Even  if 
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A  :  g  t-+  3j@,  Pp  depends  on  the  coordinate  system  chosen  in  5td,  i.e.,  it  fails  to 
have  the  desirable  property  that 


PpTg  =  TPfq 

for  all  g  6  II  and  all  ‘reasonable’  changes  of  variables,  T. 

7.  Example 

Suppose  that  A  :  g  •->  <;|e  with  0  =  {si#, . . . ,  sn0}  for  some  nontrivial  9  and 
some  (scalar)  n-set  {si,...,s„}.  Then,  for  each  k, 

A  Xk  =  [sl...,stt[6°:\a\  =  k). 

This  shows  that  each  AA*  has  rank  1.  Therefore,  regardless  of  how  the  segmented 
row-echelon  form  R  =  ( Rl} )  for  A[X0, . . .]  is  obtained,  each  has  just  one 
row,  and  this  row  is  a  scalar  multiple  of  the  vector  (9a  :  |a|  =  k).  Consequently, 
Vk  =  A'fcWfcfc  has  just  one  column,  namely  the  polynomial 

Pk  :=  Y. 

M=k 

Tliis  is  unfortunate,  for  the  following  reason.  If  6  =  q  for  some  i,  then  F  consists 
of  all  polynomials  of  degree  <  n  in  the  ith  indeterminate  only,  and  that  is  good. 
However,  if  9  has  more  than  one  nonzero  entry,  then  F  usually  differs  from  the 
natural  choice  in  that  case,  namely  the  space  of  all  polynomials  of  degree  <  n 
which  are  constant  in  any  direction  perpendicular  to  9.  For  F  to  be  that  space, 
we  would  need  to  choose  each  pk  as  some  scalar  multiple  of 

W=fc 

e.g.,  as  the  polynomial 

X  t-4  ( 9-x)k . 

This  means  that,  instead  of  W^k  =  R'kk'  we  should  choose 

Wkk  =  n* 1  J?'fcfc,  (7.1) 

with 

nfc  :=  diag(a! :  |a|  =  k).  (7.2) 

The  space  F  resulting  from  this  choice  is  the  least  choice,  i.e.,  the  space 

nA  =  0 ran J fcll*1  =:  ran[F0, VltV3,...]  (7.3) 

k 
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of  de  Boor  and  Ron  [8].  In  particular,  it  is  the  space  11©  of  [5]  in  case  A  =  |©. 
Correspondingly,  we  denote  by 


Pa,  resp.  P© 

the  resulting  interpolation  projectors. 

To  be  sure,  [8]  arrived  at  this  choice  in  a  completely  different  way,  as  the 
linear  span  of  ail  the  least  terms  of  elements  of  a  certain  space  of  formal  power 
series.  To  make  the  connection  between  (7.3)  and  the  definition  of  IIa  in  [8],  let  vt 
be  one  of  the  columns  of  [Vo,  Vj,  Vjj, . . .].  Then  vt  is  of  the  form 

v,=  ^(r/a!R(i,a) 

M=fc 


for  some  k  and,  for  that  k,  R(i,a)  =  0  for  all  )a|  <  k.  This  implies  that  vt  =  /ti 
(see  the  Introduction  for  the  definition  of  f±),  with  /,  the  formal  power  series 

/i:=^()>!f?(i,a). 

a 

Further,  R  =  M~l  A  A',  for  some  (invertible)  M.  Consequently,  the  /,  form  a  basis 
for  the  space  G  spanned  by  the  formal  power  series 

j  =  i, . . . ,n, 

a 

with  A  j  the  jth  row  of  the  data  map  A.  Further,  if  g  is  any  element  of  G  with 
deg<7;  =  fc,  then,  by  Lemma  5.1(ii),  54,  is  necessarily  in  the  range  of  Xk&kl  R'kk. 
So,  altogether,  this  shows  that  the  right-hand  side  of  (7.3)  is,  indeed,  the  linear 
span  of  all  the  leasts  of  elements  of  G,  hence  equal  to  IIa  as  defined  in  [8]. 

Taking,  in  particular,  A  =  |e,  hence  A j  =  5g  (i.e.,  point  evaluation  at  9)  for 
some  9  6  0,  we  get 


^  =  £OQ/a!0Q  =  £(-0)°/a!  =  Cfl> 
a  a 


thus  showing  that,  for  this  case,  (7.3)  does,  indeed,  give  the  space  spanned  by  all 
the  leasts  of  linear  combinations  of  exponentials  eg  with  9  €  0. 

It  has  been  pointed  out  (after  a  talk  on  this  material,  by  someone  whose 
name  I  never  learned)  that  the  simple  choice  (6.6)  actually  provides  (7.1)  if  Xk 
were  changed  to 


:=[*-►  Y[x{p{i)):p  6  <*}*]> 


i.e.,  if  the  monomials  were  treated  as  if  the  indeterminates  did  not  commute,  with 
a  corresponding  change  of  fi*,  to  the  scalar  matrix  k!. 


8.  The  least  choice 

The  least  choice,  (7.3),  not  only  provides  the  ‘right’  space  in  case  ©  lies  on 
a  straight  line,  it  has  a  rather  impressive  list  of  properties  which  are  detailed  in 
[5,7,8].  Of  these,  the  following  is  perhaps  the  most  unusual.  Its  statement  uses 
the  notation  pf,  for  the  leading  term  of  the  polynomial  p,  as  introduced  in  the 
Introduction.  Further,  if  p  =  Xc,  then 

p(D)  :=  ^D“c(q), 


with 


Da  =  D“(U 


"ud  > 


and  D{  the  derivative  with  respect  to  the  ith  argument. 

Result  8.1  ([7]). 

ne  =  P)  ker Pt(D). 
p|e=° 

A  proof  is  provided  in  [8],  as  part  of  a  more  general  argument.  A  direct  proof 
can  be  found  in  [4]. 

It  is  informative  to  consider  this  result  in  the  light  of  elimination,  as  a  further 
motivation  for  the  choice  (7.1)  for  the  Wkk-  For  this,  we  continue  to  denote  by 
MR  any  factorization  of  A[A'0,  Xi, . . .]  into  an  invertible  M  and  a  correspondingly 
segmented  row-echelon  form. 

Lemma  8.2.  Let  p  €  Tlk\n<k  with  pT  =:  Xkc.  Then,  there  exists  q  6  II<*  with 
A (p  -  q)  —  0  if  and  only  if  c  1  ran  R'kk . 

Proof.  Suppose  that  cl  ran  R'kk.  Then  Rkkc  =  0,  hence,  by  (iii)  of  Lemma 
5.1,  there  exists  b  so  that  (6,  c)  €  ker R<k-  This  implies  that  (6,c)  6  ker  AX<jt,  i.e., 
p  :=  X<k  (h,  c)  is  in  ker  A. 

Conversely,  if  p  is  of  (exact)  degree  k  and  in  ker  A,  then  p  =  X<k  (6,  c)  for 
some  b,  c  with  pf  —  Xkc,  and  0  =  Ap  =  hX<k  ( b,c )  =  MR<k  { b,c ),  hence,  since 
M  is  invertible,  also  R<k  (b,c)  =  0,  in  particular  Rkkc  =  0,  i.e.,  cl  ran  R'kk.  □ 

This  lemma  characterizes  ran  R'kk  as  the  orthogonal  complement  of  the  set 

{c  :  3{p  €  ker  A}  pt  =  Xfcc} 

of  leading  coefficients  of  polynomials  of  exact  degree  k  in  the  kernel  of  A.  Since 

iIa,*  :=  nA  n  n£  =  xkfiklmkk, 

this  provides  the  following  characterization  of  the  fcth  homogen.  is  component  of 
nA,  i.e.,  of  n^fc. 


15 


Proposition  8.3. 

n£,fc  =  1$  fl  f|{kerpt(Z>)  :  p  €  Uk\Il<k,  A p  =  0}. 

Proof.  Assume  that  q  is  a  homogeneous  polynomial  of  degree  k,  i.e., 
q  =:  Xkc,  and  consider  the  corresponding  differential  operator  q(D).  Then  TI<k  C 
ker q(D)  trivially.  Hence,  for  any  g  6  n*\n<jt,  setting  =:  Xka,  we  have 

q(D)g  =  q(D)(gr)  =  c(a)a\a(a). 

\a\=k 


In  particular,  with  Clk  as  given  in  (7.2), 

q(D)g  =  0  <=>  c±nfca. 

If  now  g  is  an  element  of  II, \,  then  =  Xk  lRkkb  for  some  b,  therefore 

Clka  —  R'kkb- 

Consequently,  IIa  n  lit  C  ker  q(D)  if  and  only  if  c±  ran  R'kk,  i.e.,  if  and  only 
if  Rkkc  =  0,  i.e.,  by  Lemma  8.2,  if  and  only  if  q  is  the  leading  term  of  some 
polynomial  in  ker  A.  □ 

Note  that,  for  the  case  A  =  |e,  Result  8.1  is  much  stronger  than  Proposition 
8.3.  For,  in  this  special  case,  not  only  is  IIq  k  annihilated  by  Pj(D)  in  case  p  6 
nfcXn^j  and  p|e  =  0,  but  all  of  He  is  annihilated  by  such  p^{D).  In  particular, 

Pt(D)XjiiJlR'jj  =  0 

for  any  j. 

The  reason  for  this  much  stronger  result  in  the  case  A  =  je  is  that,  in  this 
case,  ker  A  is  a  polynomial  ideal  and,  consequently,  IIa  Is  £Mnvariant,  i.e.,  /  6  IIa 
and  a  6  ZZ^.  implies  Daf.  See  [8]  for  details. 


9.  Construction  of  an  H-basis  for  a  polynomial  ideal  with  finite  variety 

Customarily,  a  polynomial  ideal  J  is  specified  by  some  finite  generating  set 
for  it,  i.e.,  by  describing  it  as 

i  =  £n9 

g€G 

for  some  finite  G  C  II.  Given  this  description,  it  is  nontrivial  to  determine  whether 
or  not  a  given  /  €  n  belongs  to  I,  except  when  G  is  an  H-basis  for  I,  meaning 
that,  for  every  fc, 

xk  inn*  =  y  ’ nfc-deggg. 

g€G 
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Such  a  generating  set  is  also  called  a  Macaulay  basis  or  a  canonical  basis  for 

I. 

Indeed,  if  G  is  an  H-basis  for  X,  then,  with 

1°  :=  {pt :  P  6  1} 

the  associated  homogeneous  ideal,  we  have 

^:=I°nn2  =  £n°_degsffT,  all  k. 

g€G 


Thus,  /  €  J  if  and  only  if 

ft  =  ]£  Pg9t 

g€G 

for  some  pg  €  n°_deg  9,  and  in  addition,  for  these  pg, 

f  Pg9  €  T- 

geG 

Since  the  last  condition  involves  a  polynomial  of  degree  <  deg  /,  this  leads  to  a 
terminating  recursive  check. 

If  an  ideal  is  not  given  in  terms  of  an  H-basis,  then  it  is,  in  general,  nontrivial 
to  construct  an  H-basis,  except  when  the  direct  summands  TPk  are  known  in  the 
sense  that  bases  are  known  for  them. 

In  that  situation,  one  can  construct  an  H-basis  G  for  I  by  constructing 

Gk  :=Gnnt 

inductively,  i.e.,  for  k  =  0, 1, 2, . . .,  so  that 

Tfc  —  ^  ^  Il/i-degjffi  (9-1) 

seGk 


as  follows. 

Assume  that  we  already  have  Gk- i  in  hand  (as  is  surely  the  case  initially  for 
any  k  with  Z£_j  =  {0},  e.g.,  for  k  =  0).  We  claim  that  the  choice 

Gk  Gk- 1  UB, 


with  Belt  so  that  [Bf]  is  a  basis  for  a  linear  subspace  in  Z£  complementary  to 

^fc-degjS't* 

seGt_i 
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does  the  job,  i.e.,  satisfies  (9.1). 

Indeed,  if  p  €  Ik,  then  €  Z$,  hence 

Pi  =  [Btl“  +  ps9l 

J6G4_, 


for  some  a  €  CB  and  some  pg  €  n°_degJ,  all  g  €  Gk-i-  However,  then 
p-[B]a-  P99  €  Ik- 1=  nfc-i-degs3> 

0€Gi-i  seGt-! 

therefore 

p  €  Y,  n0$  +  £  nfc  -degjSi 

S6B  9€Gt., 

which  was  to  be  shown. 

Note  that  the  H-basis  constructed  is  minimal  in  the  sense  that  any  proper 
subset  of  it  would  generate  a  proper  subideal  of  2. 

An  appropriate  B  can  be  constructed  by:  (i)  starting  with  C  for  which  [Cf] 
is  a  basis  for  Zj?_j  (as  constructed  in  the  preceding  step),  (ii)  selecting  from  the 
columns  of  [(),Jc  :  c  6  C,j  =  1  a  column  map  C  which  is  maximal  with 

respect  to  having  [CT]  1-1,  (iii)  extending  this  to  a  column  map  [C,  B]  with  (B  Cl 
and)  ,  i?f]  a  basis  for  Z£. 

Of  course,  this  requires  that  one  have  in  hand,  for  each  fc,  aBt  Cl  for  which 
[Bkx]  is  a  basis  for  Zj),  and  it  is  in  this  sense  that  one  needs  to  “know”  the  l£. 
This  information  is  easy  to  derive  in  case  I  is  an  ideal  with  finite  codimension  or, 
what  is  the  same,  with  finite  variety. 

It  is  well-known  (see,  e.g.,  de  Boor  and  Ron  [6]  for  a  detailed  retelling  of  the 
relevant  facts)  that  a  polynomial  ideal  has  finite  codimension  exactly  when  it  is  of 
the  form  ker  A,  with  A'  =:  [Ai, . . . ,  An]  a  basis  for  a  subspace  of  II'  of  the  form 

£  5ePe{D)  =  {£  5ePg{D)  :peePg}  (9.2) 

see  see 

with  0,  the  variety  of  the  ideal,  a  finite  subset  of  Cd,  Sg  :  p  t-t  p(0)  the  linear 
functional  of  evaluation  at  0,  and  each  Pg  a  ^-invariant  finite-dimensional  linear 
subspace  of  II.  To  put  it  differently,  a  correct  polynomial  interpolation  scheme 
with  data  map  A  €  L(TI,  lRn)  is  ideal  (in  the  sense  of  Birkhoff  [1])  if  and  only  if  it 
is  Hermite  interpolation,  i.e.,  its  set  of  interpolation  conditions  is  of  the  form 
(9.2). 

Of  course,  if  1  =  ker  A  for  some  explicitly  known  map  A  €  L(II,IRn),  then 
membership  in  1  of  a  given  /  €  II  is  trivially  testable:  simply  compute  A/.  On  the 
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other  hand,  an  H-basis  for  such  an  ideal  is  likely  to  have  its  uses  in  the  construction 
of  error  formulae  for  the  associated  polynomial  interpolation  schemes. 

In  any  case,  if  I  =  ker  A  for  some  onto  A  €  £,(11,  IR"),  then 

l=Xke*R, 


with 

AX  =  MR 

any  particular  factorization  of  AX  with  a  1-1  M. 

In  particular,  let  R  =  be  a  segmented  row-echelon  form,  segmented 
corresponding  to  the  segmentation  X  =  (A'0,Xi, . . .].  By  Lemma  8.2,  p  6  11“  is 
in  iff  c  :=  Xklp  €  kerf?*./..  Consequently,  any  basis  Kk  for  ker Rkk  provides 
a  basis,  A I* A'*,  for  l£.  Since  R  =  (RtJ)  is  in  segmented  row-echelon  form,  each 
Rjj  is  onto  while  RtJ  =  0  for  i  >  j;  therefore  we  can  find  Kq,  . . . ,  Kk- 1  so  that 
Hj<fc  RvKi  ~  —RtkKk  for  all  i  <  k,  while  this  holds  trivially  for  all  i  >  k.  This 
implies  that  Bk  :=  YL}<k  XjKj  is  a  column  map  into  1,  with  B ^  a  basis  for  2^. 

If  the  rows  of  R^k  are,  as  in  [7],  constructed  as  an  orthonormal  basis  (of  the 
row  space  of  a  certain  matrix  and  with  respect  to  some  inner  product),  then  a 
convenient  choice  for  Kk  is  the  completion  [R'kk,  Kk\  of  R'kk  to  an  orthonormal 
basis  for  all  of 


10.  The  construction  of  rules 

The  least  rule  for  fi  €  II'  from  A  €  L(II,  1R7*)  is,  by  definition,  the  linear 
functional 

»P\- 

Since  P\  —  V{i\V)~l  A  (with  V  some  basis  for  Ha ;  see  Lemma  3.2),  we  can  also 
write  the  least  rule  explicitly  in  terms  of  the  rows  Aj  of  A  as 

pP\  =  a  *  A  =  ^2  a(i)Aj,  (10.1) 


with 

a  :=  jiV(AV)-1. 

Since  our  basis  V  =  [Vo,  V), .. .]  =:  XW  for  II a  has  the  simple  segments 
14  =  XknklR'kk,  k  =  0,1,2,..., 

(see  (7.3))  it  is  not  hard  to  apply  the  inverse  of  the  matrix  ( AV )'  =  (AXW)'  = 
(MRW)'  to  the  vector  pXW.  In  fact,  the  resulting  work  can  be  carried  out  right 
along  with  the  construction  of  the  factorization  if  one  is  willing  to  record  M-1  in 
some  convenient  form,  as  we  now  show. 
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The  idea  is  to  apply  the  program  described  in  [7],  not  just  to  the  matrix 
AX  but  to  the  ‘augmented’  matrix  [A;  p]X.  (The  choice  of  the  word  ‘augmented’ 
here  is  quite  deliberate;  it  stresses  the  fact  that,  what  is  about  to  be  described,  is 
nothing  than  the  dual  of  the  standard  procedure  of  applying  Gauss  elimination  to 
the  augmented  matrix  [A,  6]  in  order  to  compute  a  solution  of  A?  =  6.)  For  this 
to  work,  this  additional,  last  row  must  never  be  used  as  pivot  row.  To  recall,  the 
program  in  [7]  constructs  each  R'kk  as  an  orthogonal  basis  for  the  row  space  of  a 
certain  working  array,  call  it  Bk,  whose  rows  consist  of  the  X^-part  of  each  row 
not  yet  used  as  pivot  row.  In  other  words,  the  columns  of  Bk,  i.e.,  the  entries  in 
each  row  of  Bk,  are  naturally  indexed  by  a  with  |a|  =  k.  Orthogonality  is  with 
respect  to  the  inner  product 


(b,c)k  :=  btlklc. 


(10.2) 


The  construction  consists  of  making  all  rows  in  Bk  not  yet  used  as  pivot  rows 
orthogonal  to  those  being  used  as  pivot  rows.  Without  that  additional  last  row, 
work  on  Bk  ceases  once  all  rows  not  yet  used  as  pivot  rows  are  zero.  With  that 
additional  row,  work  on  this  segment  still  ceases  when  all  rows  not  used  as  pivot 
rows  are  zero,  except  for  that  additional  row.  This  may  leave  that  additional  row 
nonzero,  but  it  is  certain  to  leave  it  orthogonal  to  the  rows  of  the  resulting  Rkk 
with  respect  to  the  inner  product  (10.2). 

At  the  end,  we  obtain  the  factorization 


A' 

’  M  O' 

R 

X  = 

m 

u 

with  AX  =  MR  the  earlier  factorization,  and  with  u  the  segmented  sequence 
(u0,  tii, . . .)  such  that 

Uk-nklR'kk  =  0 

for  every  nonempty  Rkk ■  One  verifies  that 


M  O' 

-l 

’  M-1 

O' 

TO  1 

1 

(Here,  m  is  a  1  x  (m  -  l)-matrix  rather  than  a  sequence,  hence  I  write  mM  1 
rather  than  the  (undefined)  This  implies  that 

(-mAT’A  +  n)V  =  (— mA/-1A  +  p)XW  =  £  ukQklR'kk  =  0. 

k 

Consequently,  mM~l AV  =  pV,  hence  the  desired  coefficient  vector  a  for  (10.1) 
equals  the  one  row  of  mM-1. 
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11.  A  second  look  at  Gauss  elimination  by  segments 


Here  is  a  second  look  at  Gauss  elimination  by  segments,  from  the  point  of  view 
of  the  data  map  A  6  L(II,IRn).  This  discussion  is  applicable  to  any  (segmented) 
matrix  A  =  [Ao,  Ai, . . .]  since  an  arbitrary  matrix  A  €JFnxJ  can  be  thought  of  as 
such  a  product  AX  of  a  row  map  A  €  L(S,  Fn)  with  a  column  map  X  G  L(JFJ ,  5) 
in  at  least  two  ways: 

(1)  S  =  Fn  and  A  =  id,  hence  X  =  A  =  [A(:,  j)  :  j  G  J). 

(2)  5  =  F7  and  V  —  id,  hence  A  =  A  G  L(5,  Fn). 

At  the  same  time,  it  is  to  be  hoped  that  elimination  by  segments  will  be  found 
of  help  in  other  applications  where  it  is  more  natural  to  do  elimination  by  cer¬ 
tain  segments  rather  than  by  columns.  For  example,  a  process  requiring  a  certain 
amount  of  column  pivoting  might  be  a  good  prospect. 

For  this  reason,  in  this  section,  X  =  [A'o,  At , . . .]  is  any  segmented  column 
map,  from  IF^  :=  (J  -*  F)o  into  some  linear  space  S,  with  J  some  set,  and  A  is 
any  onto  data  map  or  row  map  on  5,  and  the  matrix  of  interest  is 

A  :=  AX  =  [AA'o,  A.Yi , . . .]  =:  [Ao, 


as  before. 

We  have  already  made  use  of  the  following  facts:  With  At, . . . ,  A„  the  rows 
of  the  data  map  A  €  L(5,Fn),  the  dual  map  A'  €  L(Fn,S')  for  A  is  of  the  form 

A'  =  (At, ...  ,AnJ :  F"  ->  S' :  a  £a(i)A,. 

i 

Since  A  is  onto  by  assumption,  A'  is  1-1,  hence  a  basis  for  its  range. 

For  each  k,  let  Ak  be  a  data  map  whose  rows,  when  restricted  to  Sk  := 
ran  A'<fc,  provide  a  basis  for  the  linear  space 

Lk  :=  {A|Sk  :  A  e  L,  AX<fc  =  0}. 

Equivalently,  Ak  is  any  onto  data  map  with  its  dim  Lk  rows  taken  from  ran  A',  for 
which  A kXKk  =  0  while  A is  onto. 

However  this  is  actually  done  numerically,  we  will  have  constructed  a  basis 

[Ai-  :  fc  =  1, 2, . . .] 

for  ran  A'.  As  A'  is  1-1,  hence  a  basis  for  its  range,  this  means  that  we  have,  in 
effect,  constructed  an  invertible  matrix  M1  €  lRmxm  so  that 

A'  =  K:*  =  0,1,2,...]M'. 


Thus, 


A X  =  A-  MR, 
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with 


R  :=  [Ao;  A*;  A2;  - .  -1  [Xo,  X*,  X2,  •  • .]  ='.  (R,j  :  i,  j  —  0, 1 ,2, . . .) 
a  block  upper  triangular  matrix  since 

Rij  =  A,X, 

is  trivial  tor  i  >  j. 

Consider  now,  in  particular,  the  diagonal  blocks, 


Rkk  =  AjfcXfc. 


If  A  is  any  linear  functional  in  ran  A'  which  vanishes  on  II<fc,  then  there  is  exactly 
one  coefficient  sequence  c  so  that  A  =  A'fcc  on  11*...  In  other  words,  the  rows  of  /?*.*. 
provide  a  basis  for  the  linear  subspace  {(A()a  :  |a|  =  k)  :  A  G  ranA\AX<fc  =  0}. 
In  particular,  R'kk  is  1-1,  and  its  range  is  independent  of  the  particular  choice  of 
the  map  A'fc.  This  recovers  (i)  of  Corollary  5.2  in  this  more  general  setting. 
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